Entropy driven atomic motion in laser-excited bismuth 
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We introduce a thermodynamical model based on the two-temperature approach in order to fully 
understand the dynamics of the coherent Ai s phonon in laser-excited bismuth. Using this model, 
we simulate the time evolution of (111) Bragg peak intensities measured by Fritz et al [Science 
315, 633 (2007)] in femtosecond X-ray diffraction experiments performed on a bismuth film for 
different laser fluences. The agreement between theoretical and experimental results is striking not 
only because we use fluences very close to the experimental ones but also because most of the model 
parameters are obtained from ab-initio calculations performed for different electron temperatures. 

PACS numbers: 61.80.Ba, 63.20.dk, 71.15.Pd, 78.47.J- 



The development of optical pump-probe experiments 
has revolutionised physics by allowing the study of 
non-equilibrium processes on a femtosecond time scale. 
Indeed, an ultrashort laser pulse can cause dramatic 
changes in the physical properties of covalently bonded 
solids and induce ultrafast melting, phase transitions 
or coherent phonons [l|. Coherent excitation of optical 
phonons is a general phenomenon which has been ex- 
tensively investigated in bismuth where the reflectivity 
change following electronic excitations exhibits damped 
oscillations arising from atomic motions corresponding to 
the Ai g mode [2]. Unfortunately, it is not possible to di- 
rectly relate the amplitude of the reflectivity signal to the 
atomic positions. The goal of measuring atomic position 
as a function of time remained elusive until the break- 
through of time resolved X-ray diffraction experiments 
which allow atomic motions to be followed stroboscop- 
ically with picometer spatial and femtosecond temporal 
resolution [3]. Different femtosecond X-ray diffraction ex- 
periments have been conducted to probe the dynamics of 
the Ai g phonon in laser-excited Bi043- I n this letter, we 
try to understand and to reproduce the time evolution of 
(111) Bragg peak intensities measured by Fritz et ai||. 

Bismuth is a semi-metallic solid, whose structure can 
be derived from the face-centered-cubic lattice by a weak 
rhombohedral distorsion of the cubic unit cell. The crys- 
tal basis consists of two atoms located along the body 
diagonal of length c at ±uc, where u is the Ax g phonon 
coordinate. At 300 K, the equilibrium value of u is 
u eq ~ 0.234j8|. It is now well established that the sudden 
change of u eq following laser excitation explains the gen- 
eration of coherent phonons in Bi. The time evolution of 
u can be obtained from the displacive excitation of co- 
herent phonons (DECP) modelQ. This model, based on 
the idea that u eq increases linearly as a function of the 
excited carriers density n and that n decreases exponen- 
tially as a function of time, however does not capture all 
aspects of the phonon and electron dynamics. Two ap- 
proaches based on ab-initio calculations [9l [lfl| have been 
proposed to describe the non-equilibrium electron distri- 
bution. Murray et al^ assumed that electrons and holes 
can be described by two Fermi-Dirac (FD) distributions 



with the same temperature but different chemical poten- 
tials, and obtained a good agreement with experiments 
concernin g th e early oscillations of the Ai ff phonon[5j. Zi- 
jlstra et al[jj^ studied the coupling between the Ai s and 
E s modes and assumed that electrons and holes can be 
described by a unique FD distribution. They claimed 
that a two-chemical potential model does not provide 
a realistic description of the electron dynamics in Bi, 
at least for high fluences. While such a claim seems 
to be supported by recent work based on a two-fluid 
model [6j, t he relevance of the two approaches remains 
unclearjjlOJ]. 

In this work, we consider a model based on the two- 
temperature approach [u| where the electron system is 
locally described by a single FD distribution, the A lg 
phonon mode obeys a classical equation of motion and 
all the remaining modes are lumped in a unique phonon 
bath with specific heat C/. We assume that both the 
electron temperature T e and lattice temperature T\ are 
well defined at each time because of electron-electron and 
phonon- phonon interactions. Making a local energy bal- 
ance, the heat respectively received by the electron and 
phonon systems is given by: 

dS dU fedu d ( dT e \ 

(1) 

. 4Mc 2 / du\ 2 

-GoK-vn — f-) (2) 



Ci 



m 
dt 



where P is the energy deposited by the laser pulse in 
the electron system per unit volume and time, S and 
U are electronic entropy and energy per unit volume, / 
is the force acting on the Ai g phonon, n e is the electron 
thermal conductivity, Go is the electron-phonon coupling 
constant which governs the heat transfer from the elec- 
tron system to the lattice, v is the unit cell volume, M 
is the mass of a Bi atom and r is the lifetime of the Ai g 
phonon. Since the laser spot diameter is much larger 
than the penetration depth we only consider the spatial 
dependence of the parameters in the direction perpendic- 
ular to the surface, labelled as z. One should notice that 
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all parameters of Eq. [T] and [2] depend on T e and u, both 
of which depend on the depth z. These equations have to 
be solved together with the phonon equation of motion : 
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where the last term describes the energy transfer from 
the coherent phonon to the incoherent phonon bath, and 
where the force / is defined by: 



v dU 
c du 
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In order to evaluate the force / acting on the coher- 
ent phonon and other quantities involved in our model, 
we performed ab-initio calculations for different elec- 
tron temperatures T e and phonon coordinates u. All 
these calculations were done within the framework of 
the local density approximation (LDA) for the exchange- 
correlation functional to density functional theory (DFT) 
using the ABINIT code 14] . Spin-orbit coupling was 
included and an energy cut-off of 15 Hartree in the 
planewave expansion of wavefunctions together with a 
12x12x12 kpoint grid for the Brillouin zone integration 
were used. These parameters ensure the convergence of 
both vibrational and electronic properties 15j. The ex- 
perimental lattice parameters at room temperature are 
used[|| and volume changes due to an increase of the 
lattice temperature 7} are neglected. 

First, we calculated both energy U and entropy S as 
a function of u and T e and hence obtained U(S,u). Fig- 
ure [T] a shows that the equilibrium position u eq increases 
from 0.233 to 0.25 as S increases from to ~ 2.0 kg per 
unit cell. This finding confirms the displacive mechanism 
of phonon generation and is qualitatively in line with the 
behaviour observed by Murray et a/[9j] usin g a two chem- 
ical potential approach or by Zijlstra et al 10] using the 
same type of approach as in the work reported here. If 
we assume that the laser pulse energy is deposited homo- 
geneously at the end of the pulse (n e — > oo), Eq. [ljshows 
that the electronic entropy S becomes a constant of mo- 
tion provided that the energy transfer from the electron 
system to the lattice is frozen (Go = 0). The solution 
of Eq. [3] when t — > oo (undamped dynamics) then pro- 
vides the phonon frequency v of the Ai s mode for a given 
entropy S which depends on the laser pulse fluence. Fig- 
ure [JJb shows that v decreases from 2.93 Thz to Thz 
when S increases from to ~ 1.37 kg per unit cell. This 
red shift is still in qualitative agreement with previous 
theoretical studies (9l. [Toj. 

We have also calculated the lattice specific heat Ci 
within the harmonic approximation [l6j]. As shown in Fig. 
[TJc, the agreement between our calculated lattice specific 
heat and the experimental datafjjj is good up to 7} ~ 300 
K. The discrepancies between theoretical and experimen- 
tal results become larger at higher temperature because 
of anharmonic interactions neglected in our calculation. 
In our simulations, the initial lattice temperature T was 
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FIG. 1: (color online) a/ u eq as a function of electronic en- 
tropy S. b/ Ai 9 phonon frequency v (in Thz) as a function 
of S. The star indicates the Raman frequency measured at 
300 K[l7|. c/ The calculated constant- volume lattice spe- 
cific heat Ci (solid line) of Bi compared to experimental data 
(open circles) from Ref. [17| for lattice temperatures Ti up to 
the melting temperature Tf = 544 K. d/ Calculated electron 
specific heat C e as a function of electron temperature T e for 
w=0.233 (black circles) and -u=0.24 (violet triangles). 



set to 300 K. As T° > 6 D where 9 D = 119 K is the Debye 
temperature, the temperature dependence of Ci can be 
neglected. Indeed, C;(T°) is only 1 % smaller than the 
value given by the Dulong-Petit law shown as a dashed 
line in Fig. [TJc. The electron temperature dependence of 
Ci can also be safely neglected. 

Finally, we calculated the electron specific heat of Bi, 
C e = dU/dT e , as a function of electron temperature 
T e for different values of u. As expected, the temper- 
ature dependence of C e is linear at low electron tem- 
peratures. Figure [TJd shows that a significant posi- 
tive deviation from the linear behaviour occurs when 
T e becomes larger than ~ 300 K and that this devia- 
tion increases when u increases, i.e. when Bi goes to- 
ward the hypothetical metallic state corresponding to 
u = 0.25. The knowledge of C e (T e ) is crucial, not 
only because it gives a correct estimate of the rise of 
temperature in the electron system, but also because it 
is related to the electron thermal conductivity n e . In- 
deed, assuming that the total electron scattering rate is 
dominated by the electron-phonon scattering rate, the 
electron thermal conductivity can be approximated by 
K e (T e ,T t ) = Ko x (C e (T e )/C e (T°j) x (T°/T,)|ll El 
where the dependence on u has been omitted and where 
kq = 11 W.m^.K -1 is the experimental value at room 
temperature [2fJ| . 

We can now simulate the time-resolved X-ray diffrac- 
tion experiments performed by Fritz et al^ on a Bi film 
of thickness L — 50 nm excited by a near-infrared laser 
pulse whose full width at half maximum (FWHM) is 
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FIG. 2: (color online) Normalized intensities of the (111) 
Bragg peak measured by Fritz et a/[f| as a function of time 
delay t between the laser pulse and X-ray probe for differ- 
ent absorbed fluences (symbols) compared to the theoreti- 
cal calculation (solid lines). The measured fluences are 0.7 
(black circles), 1.2 (red triangles), 1.7 (blue squares) and 
2.3 mj.cm -2 (magenta losanges). The parameters used in 
our simulations are given from top to bottom for curves la- 
belled 1 to 4: the absorbed fluences F a bs are 0.64 (black), 
1.26 (red), 2.18 (blue) and 2.47 mj.cm -2 (magenta); The ef- 
fective electron-phonon coupling constants Go are 0.84, 1.15, 
1.40 and 1.45xl0 16 W.m^.K" 1 and the decay times r are 
0.87, 0.64, 0.44 and 0.25 ps respectively. The uncertaintitics 
for r are quite large for curves 1 and 2 because of the limited 
duration of the experiments. The fitted value of the FWHM 
of the X-ray pulses is tx = 205 fs. 
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FIG. 3: (color online) a/u as a function of time delay t for the 
parameters given in the caption of Fig. [2] b/Electron tem- 
perature T e as a function of time delay t. The inset displays 
the variation of the lattice temperature Ti as a function of t. 
The curves labelled 1 to 4 are related to the curves with the 
same label in Fig. [2] 
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70 fs. The source term P(z,t) in Eq. [T] is given 
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where l p = 14 nm is the penetration depth at wavelength 
of 800 nm[T7j and F a b s is the absorbed fluence. The so- 
lution of the coupled differential equations [TJ [2] and [3] for 
a set of the two unknown physical parameters Go and r, 
which are assumed to remain constant on the time scale 
of the experiment for a given fluence, provides u(z, t) and 
T e (z, t). As will be seen later, the spatial dependence of 
U and T e can be neglected only ~ 100 fs after the ar- 
rival of the laser pulse on the Bi surface. Therefore, the 
normalized intensity of the (111) Bragg peak is given by 
I(t)/I(0) = cos 2 [67ra(£)]/ cos 2 [67ru(0)] and is convoluted 
with a Gaussian with FWHM tx to account for the tem- 
poral resolution of the experiment {BJ. The results of a 
least-squares fit of our model to the experimental data 
are shown as solid curves in Fig. [2] for four absorbed 
fluences. The parameters F a b Sl Go, t and tx are given 
in the caption. The absorbed fluences used in our calcu- 
lations are very close to the experimental fluences with 
the largest deviation occuring for the curve labelled 3 
in Fig. f2] The electron-phonon coupling constant Go 
increases from 0.84 to 1.45 xlO 16 W.m~ 3 .K _1 and the 
decay time r of the coherent phonon decreases from 0.87 



to 0.25 ps as the theoretical fluence increases from 0.64 
to 2.47 mJ.cm~ 2 . The reported electron-phonon cou- 
pling constants are small compared to usual values in 
metals 18] because the density of states at the Fermi level 
is very small for Bi. Figure [2] shows that the oscillations 
of the (111) Bragg intensities are nicely reproduced by 
our quasi-isentropic model with the exception of curve 
3, where F a b s is overestimated with respect to the ex- 
perimental fluence. It is worth remarking that curves 
3 and 4 are very close to each other while the experi- 
mental fluences differ by 0.6 mJ.cm" 2 . Interestingly, an 
upward shift of curve 3 by about 0.05 leads to a bet- 
ter agreement with our theoretical intensity for a fluence 
f a &s=l-79 mJ.cm~ 2 very close to the experimental one. 

One can now discuss the physical processes playing a 
role in the generation of coherent phonons. The laser 
energy is initially deposited in the electron system since 
the laser pulse duration is much shorter than the charac- 
teristic time for the electron-lattice energy exchange. At 
the same time, the heat diffuses in the electron system. 
The diffusivity x is approximately constant and given by 
Ko/C e (T°) ~ 100 cm 2 .s -1 . Therefore, the usual electron 
heat transport occurs and the time needed to uniformizc 
the electron temperature T e is roughly given by L 2 /\ ~ 
250 fs. Our simulations show that the spatial variations 
of T e , S and u can be safely neglected only 100 fs after 
the laser pulse maximum. The strong overheating of the 
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electron system at the end of the laser pulse is illustrated 
in Fig. (3]b. The electron temperature reaches 2231 K 
for the lowest fluence and 3724 K for the highest fluence. 
The corresponding increases of electronic entropy S are 
respectively 0.39 and 0.98 ks per unit cell. As shown 
in Fig. [T]a, the increase in entropy leads to an increase 
in u eq . The atoms thus start to oscillate around a new 
equilibrium position which evolves slowly as a function 
of time t due to heat transfer from the electron system 
to the lattice. The damped oscillations of the reduced 
coordinate u as a function of time delay t are shown in 
Fig. [31a. The initial red-shift of the phonon frequency v 
as the fluence increases is clearly visible and can be at- 
tributed to the decrease of v as S increases (see Fig. (Ub). 
Figure [3] also shows that the oscillations of the phonon 
coordinate u are accompanied by oscillations in the elec- 
tron temperature T e whose amplitude and damping grow 
with fluence. The oscillations in T e and u are in antiphase 
with respect to each other. Such a behaviour is easy to 
understand if one assumes that the electron subsystem 
undergoes an isentropic transformation (Go — > 0). As the 
metallicity of Bi is enhanced when u increases, the elec- 
tronic entropy S at constant electron temperature would 
increase up to a maximum value for the largest value of u. 
In order to compensate for the increase of S, T e decreases 
and reaches its minimum value when u reaches its maxi- 
mum value. Thus, T e oscillates at the same frequency as 
u but in antiphase. Interestingly, such a behaviour paral- 
lels that of a gas enclosed in an insulating piston chamber. 
Assuming that the gas undergoes an isentropic transfor- 
mation, the piston when displaced from it's equilibrium 
position, performs an oscillatory motion while the gas 
temperature undergoes a similar oscillatory motion, but 
in antiphase. In Bi, the electronic entropy S does not 
remain constant after the arrival of the laser pulse but 



decays slowly due to energy exchange between the elec- 
tron and lattice subsystems. The inset of Fig. (3]b shows 
the evolution of the lattice temperature T\ as a function 
of time delay t. On the time scale of the experiments 
(~ 2.5 ps), the increase in 1) ranges from 34 K for the 
lowest fluence to 110 K for the highest fluence. Although 
the increase of the lattice temperature is significant for 
the highest fluence, volume changes are expected to be 
delayed as illustrated by the shift in the Bragg angle oc- 
curing on a time scale of ~ 20 ps in experiments [5j. It's 
thus reasonable to assume that the lattice parameters 
remain constant in our simulations. 

In conclusion, we have developed a thermodynamical 
model in order to simulate the time evolution of the A\ g 
phonon coordinate following the arrival of a laser pulse 
of a given fluence on a Bi film. The intensities of the 
(111) Bragg peak measured by Fritz et ai[5j are fairly 
well reproduced by our model for fluences very close to 
the experimental ones. This success is noteworthy since 
the force acting on the coherent phonon as well as most of 
the model parameters are obtained from ab-initio calcu- 
lations. The only adjustable parameters are the effective 
electron-phonon coupling constant Go and the scattering 
time 1 /t of the coherent phonon. Our results show that 
(1) both parameters increase as the fluence increases, (2) 
the electronic heat diffusion is crucial and (3) the oscil- 
lations of the coherent phonon are accompanied by an 
oscillation in the electron temperature. From a funda- 
mental point of view, this work firmly establishes that a 
single chemical potential approach is reliable for describ- 
ing the excited electrons and provides a complete scenario 
for the generation of coherent phonons in Bi films. 
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